
pfix = vgg_get_nonhomg(X);
pfix = transpose(pfix);
%figure;
dt = DelaunayTri(pfix);
tetramesh(dt, 'FaceAlpha', 0.1);
axis equal


R = 10;
%figure;
T = delaunayn(pfix,{'Qt','Qbb','Qc','Qz'});
    % Limit circumradius of primitives
    if R < inf
        dt = TriRep(T,pfix);
        [rcc,rcc] = circumcenters(dt);
        T = T(rcc < R,:);
    end    
% Remove inner faces
F = [T(:,2:4); T(:,[1 3 4]); T(:,[1 2 4]); T(:,1:3)];
F = sort(F,2);
F = sortrows(F);
ii = find(~any(diff(F),2));
F([ii;ii+1],:) = [];

% Plot boundary faces
rotate3d on;
trisurf(F,pfix(:,1),pfix(:,2),pfix(:,3),'FaceColor','yellow','FaceAlpha',0.5);
axis equal

